# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import numpy as np
from hysop.tools.htypes import check_instance, first_not_None
from hysop.tools.misc import upper_pow2_or_3
from hysop.backend.device.codegen.base.enum_codegen import EnumCodeGenerator
from hysop.backend.device.codegen.base.struct_codegen import StructCodeGenerator
from hysop.backend.device.opencl.opencl_types import OpenClTypeGen
from hysop.constants import BoundaryCondition
BoundaryConditionEnum = EnumCodeGenerator(
BoundaryCondition, comments="Mesh boundary conditions enumeration"
)
[docs]
class MeshBaseStruct(StructCodeGenerator):
def __init__(self, typegen, vsize, typedef=True):
assert vsize in [1, 2, 3, 4, 8, 16]
name = f"{typegen.fbtype[0]}MeshBase{vsize}D"
dtype, comments = MeshBaseStruct.build_dtype(typegen, vsize)
if isinstance(typedef, bool):
if typedef is True:
typedef = name + "_s"
else:
typedef = None
super().__init__(
name=name, dtype=dtype, typegen=typegen, typedef=typedef, comments=comments
)
self.reqs["boundary_enum"] = BoundaryConditionEnum
self.vsize = vsize
[docs]
@staticmethod
def build_dtype(typegen, vsize):
tg = typegen
intn = tg.dtype_from_str(f"int{vsize}")
floatn = tg.dtype_from_str(f"fbtype{vsize}")
dtype = []
for fn in ["resolution", "compute_resolution"]:
field = (fn, intn)
dtype.append(field)
for fn in ["xmin", "xmax", "size"]:
field = (fn, floatn)
dtype.append(field)
dtype.append(("lboundary", intn))
dtype.append(("rboundary", intn))
comments = [
"Resolution of the mesh -including- ghosts",
"Resolution of the mesh -excluding- ghosts",
"Coordinates of the 'lowest' point including ghosts",
"Coordinates of the 'highest' point including ghosts",
"size = xmax - xmin",
"Left boundary conditions",
"Right boundary conditions",
]
return np.dtype(dtype), comments
[docs]
def create(
self,
name,
resolution,
compute_resolution,
boundaries,
xmin,
xmax,
size,
**kargs,
):
vsize = self.vsize
dtype = self.dtype
def extend(var, d=0):
return np.asarray(tuple(var) + (d,) * (vsize - len(var)))
tg = self.typegen
lboundary, rboundary = boundaries
lboundary = extend(lboundary, BoundaryCondition.NONE)
rboundary = extend(rboundary, BoundaryCondition.NONE)
_lboundary = [bd() for bd in lboundary]
_rboundary = [bd() for bd in rboundary]
mesh_base_vals = {
"resolution": tg.make_intn(resolution, vsize),
"compute_resolution": tg.make_intn(compute_resolution, vsize),
"lboundary": tg.make_intn(_lboundary, vsize),
"rboundary": tg.make_intn(_rboundary, vsize),
"xmin": tg.make_floatn(xmin, vsize),
"xmax": tg.make_floatn(xmax, vsize),
"size": tg.make_floatn(size, vsize),
}
var = np.empty(shape=(1,), dtype=dtype)
for k in mesh_base_vals:
var[k] = mesh_base_vals[k]
lboundary = BoundaryCondition.array_variable(
"lboundary", typegen=tg, vals=lboundary
)
rboundary = BoundaryCondition.array_variable(
"rboundary", typegen=tg, vals=rboundary
)
value = dict(
resolution=extend(resolution),
compute_resolution=extend(compute_resolution),
xmin=extend(xmin),
xmax=extend(xmax),
size=extend(size),
)
var_overrides = {"lboundary": rboundary, "rboundary": lboundary}
cg_var = self.build_codegen_variable(
name=name, value=value, var_overrides=var_overrides, **kargs
)
return (var, cg_var)
[docs]
class MeshInfoStruct(StructCodeGenerator):
def __init__(self, typegen, vsize, typedef=True, mbs_typedef=None):
assert vsize in [1, 2, 3, 4, 8, 16]
name = f"{typegen.fbtype[0]}MeshInfo{vsize}D"
bname = f"{typegen.fbtype[0]}MeshBase{vsize}D_s"
mbs_typedef = first_not_None(mbs_typedef, bname)
if isinstance(typedef, bool):
if typedef is True:
typedef = name + "_s"
else:
typedef = None
dtype, comments, ctype_overrides, reqs = MeshInfoStruct.build_dtype(
typegen, vsize, mbs_typedef=mbs_typedef
)
super().__init__(
name=name,
dtype=dtype,
typegen=typegen,
typedef=typedef,
comments=comments,
ctype_overrides=ctype_overrides,
)
for req in reqs:
self.require(req.name, req)
self.mesh_base = reqs[0]
self.vsize = vsize
[docs]
@staticmethod
def build_dtype(typegen, vsize, mbs_typedef=None):
tg = typegen
mbs_typedef = first_not_None(
mbs_typedef, f"{typegen.fbtype[0]}MeshBase{vsize}D_s"
)
intv = tg.dtype_from_str(f"int{vsize}")
floatv = tg.dtype_from_str(f"fbtype{vsize}")
i = 0
dtypes = []
def _append(dtype):
dtypes.append(dtype)
_append(("dim", np.int32))
i += 1
for fn in ["start", "stop", "ghosts"]:
_append((fn, intv))
i += 1
for fn in ["dx", "inv_dx"]:
_append((fn, floatv))
i += 1
mesh_base = MeshBaseStruct(tg, vsize, typedef=mbs_typedef)
_append(("local_mesh", mesh_base.dtype))
_append(("global_mesh", mesh_base.dtype))
i += 2
comments = [
"Dimension of the mesh",
"Index of the first local compute point in the global grid",
"Index of the last local compute point in the global grid",
"Number of ghosts in each direction",
"Space discretization",
"1/dx",
"Local mesh",
"Global mesh",
]
return dtypes, comments, None, [mesh_base]
[docs]
def create(
self, name, dim, start, stop, ghosts, dx, local_mesh, global_mesh, **kargs
):
vsize = self.vsize
if dim > vsize:
msg = f"Dim should be less or equal to {vsize}, got dim={dim}."
raise ValueError(msg)
tg = self.typegen
vsize = self.vsize
dtype = self.dtype
dx = np.asarray(dx)
mesh_info_vals = {
"dim": tg.make_intn(dim, 1),
"start": tg.make_intn(start, vsize),
"stop": tg.make_intn(stop, vsize),
"ghosts": tg.make_intn(ghosts, vsize),
"dx": tg.make_floatn(dx, vsize),
"inv_dx": tg.make_floatn(1.0 / dx, vsize),
"local_mesh": local_mesh[0],
"global_mesh": global_mesh[0],
}
def extend(var, d=0):
if isinstance(var, np.ndarray):
_var = np.full(shape=(vsize,), fill_value=d, dtype=var.dtype)
_var[: var.size] = var
return _var
else:
return np.asarray(tuple(var) + (d,) * (vsize - len(var)))
value = dict(
dim=dim,
start=extend(start),
stop=extend(stop),
ghosts=extend(ghosts),
dx=extend(dx),
inv_dx=extend(1.0 / dx),
)
local_mesh[1].name = "local_mesh"
global_mesh[1].name = "global_mesh"
var_overrides = dict(
local_mesh=local_mesh[1],
global_mesh=global_mesh[1],
)
cg_var = self.build_codegen_variable(
name=name, value=value, var_overrides=var_overrides, **kargs
)
var = np.empty(1, dtype)
for k in mesh_info_vals:
var[k] = mesh_info_vals[k]
return (var, cg_var)
[docs]
@staticmethod
def create_from_mesh(name, typegen, mesh, **kargs):
from hysop.mesh.cartesian_mesh import CartesianMeshView
check_instance(name, str)
check_instance(mesh, CartesianMeshView)
check_instance(typegen, OpenClTypeGen)
tg = typegen
dim = mesh.dim
vsize = upper_pow2_or_3(dim)
btd = f"{tg.fbtype[0]}MeshBase{vsize}D_s"
itd = f"{tg.fbtype[0]}MeshInfo{vsize}D_s"
mesh_base = MeshBaseStruct(tg, vsize, typedef=btd)
mesh_info = MeshInfoStruct(tg, vsize, typedef=itd, mbs_typedef=btd)
start = mesh.global_start[::-1]
stop = mesh.global_stop[::-1]
dx = mesh.space_step[::-1]
ghosts = mesh.ghosts[::-1]
gresolution = mesh.grid_resolution[::-1]
gcompute_resolution = mesh.grid_resolution[::-1]
gxmin = mesh.global_origin[::-1]
gxmax = mesh.global_end[::-1]
gsize = mesh.global_length[::-1]
gboundaries = mesh.global_boundaries[::-1]
lresolution = mesh.local_resolution[::-1]
lcompute_resolution = mesh.compute_resolution[::-1]
lxmin = mesh.local_origin[::-1]
lxmax = mesh.local_end[::-1]
lsize = mesh.local_length[::-1]
lboundaries = (mesh.local_boundaries[0][::-1], mesh.local_boundaries[1][::-1])
const = kargs["const"] if ("const" in kargs) else False
gmesh = mesh_base.create(
name + "_global_mesh",
gresolution,
gcompute_resolution,
gboundaries,
gxmin,
gxmax,
gsize,
const=const,
)
lmesh = mesh_base.create(
name + "_local_mesh",
lresolution,
lcompute_resolution,
lboundaries,
lxmin,
lxmax,
lsize,
const=const,
)
(var, cg_var) = mesh_info.create(
name, dim, start, stop, ghosts, dx, lmesh, gmesh, **kargs
)
return (var, cg_var)
[docs]
def build_codegen_variable(self, name, var_overrides=None, **kargs):
tg = self.typegen
if var_overrides is None:
const = kargs["const"] if ("const" in kargs) else False
var_overrides = dict(
local_mesh=self.mesh_base.build_codegen_variable(
"local_mesh", const=const
),
global_mesh=self.mesh_base.build_codegen_variable(
"global_mesh", const=const
),
)
return super().build_codegen_variable(
name=name, var_overrides=var_overrides, **kargs
)
if __name__ == "__main__":
from hysop.backend.device.codegen.base.opencl_codegen import OpenClCodeGenerator
from hysop.backend.device.codegen.base.test import _test_typegen
tg = _test_typegen("double", float_dump_mode="dec")
vsize = 4
mbs = MeshBaseStruct(tg, typedef=f"{tg.fbtype[0]}MeshBase{vsize}D_s", vsize=vsize)
mis = MeshInfoStruct(
tg,
typedef=f"{tg.fbtype[0]}MeshInfo{vsize}D_s",
mbs_typedef=mbs.typedef,
vsize=vsize,
)
cg = OpenClCodeGenerator("test_generator", tg)
# declare mesh MeshInfoStruct and its dependancies (MeshBaseStruct,MeshDirectionEnum,TranspositionStateEnum)
cg.require("mis", mis)
# create a local numpy and a codegen MeshInfoStruct variable
local_mesh = mbs.create(
"local",
(
10,
10,
10,
),
(
8,
8,
8,
),
((BoundaryCondition.PERIODIC,) * 3,) * 2,
(
0,
0,
0,
),
(
1,
1,
1,
),
(
1,
1,
1,
),
)
global_mesh = mbs.create(
"global",
(
100,
100,
100,
),
(
80,
80,
80,
),
((BoundaryCondition.PERIODIC,) * 3,) * 2,
(
0,
0,
0,
),
(
10,
10,
10,
),
(
10,
10,
10,
),
)
(np_mis, cg_mis) = mis.create(
"mesh_info",
3,
(
0,
0,
0,
),
(
1024,
1024,
1024,
),
(
1,
1,
1,
),
(0.1, 0.2, 0.3),
local_mesh,
global_mesh,
storage="__constant",
)
# declare and intialize the nested struct in the __constant address space
cg.jumpline()
cg_mis.declare(cg)
# show generated code in your editor
cg.edit()
# build generated code on every OpenCL device to check
cg.test_compile()